rm(list=ls())
gc()
require(foreign)
require(ineq)
require(dplyr)
library(Hmisc)
library(lavaan)


setwd(paste(substr(getwd(),1,regexpr("Dropbox",getwd())[1]-1),"Dropbox/coauthors/AMT_Bisbee_Larson/bisbee_larson_replication_final",sep=""))
getwd()
source("./tools/leg2.R")
###################################################################################
##
##  File Name:    protest_results.R
##
##  Input Files:  prot-simple.dta
##  Output Files: SI_prot_causal.pdf             (Figure 5)
##
##  Purpose:  This R script summarizes the causal relationship between tie dimension
##            and decision to join a protest by exploiting random assignment of 
##            elicited tie in the data. 
##
##############################################################################################

ttest.coefs <- function(beta1,beta2,se1,se2) {
  t.stat <- (beta1-beta2)/(sqrt(se1^2 + se2^2))
  pval <- 2*pt(-abs(t.stat),df = 450)
  star <- paste("(",sprintf("%.2f",round(pval,2)),ifelse(pval < 0.05,"*)",")"),sep="")
  return(star)
}


# ----------------------------------------- Causal ------------------------------------- #
dat.cause <- read.dta("./protests/prot-simple.dta")
dat.cause$prot_dim_interaction <- dat.cause$prot_dim_interaction*-1
dat.cause$prot_dim_seat_synth2_mst <- dat.cause$prot_dim_seat_synth2_mst*-1
dat.cause$prot_dim_seat_synth2_mst_rob <- dat.cause$prot_dim_seat_synth2_mst_rob*-1

dims <- c("jobsearch","contribute","garner","personalgain","grpgain","homo_pol","homo_rel","homo_educ",
          "homo_ls","pers_cris_synth","pers_succ_synth","prof_cris_synth","prof_succ_synth",
          "seat_synth2_lst","seat_synth2_mst","interaction")
coefs.biv <- ses.biv <- NA
for(i in 1:length(dims)) {
  coefs.biv[i] <- summary(lm(paste("prot_join ~ prot_dim_",dims[i]," + prot_peace",sep=""),dat.cause))$coefficients[2,1]
  ses.biv[i] <- summary(lm(paste("prot_join ~ prot_dim_",dims[i]," + prot_peace",sep=""),dat.cause))$coefficients[2,2]
}
names(coefs.biv) <- names(ses.biv) <- dims

model.cause <- paste("prot_join ~ ",paste(paste("prot_dim_",dims,sep=""),collapse=" + ")," + weakest_dum + weak_dum + strong_dum + strongest_dum + prot_peace",sep = "")
fit.cause <- sem(model.cause,data = dat.cause)
coefs.cause <- fit.cause@ParTable$est[which(grepl("^~$",fit.cause@ParTable$op))][1:16]
ses.cause <- fit.cause@ParTable$se[which(grepl("^~$",fit.cause@ParTable$op))][1:16]

ords <- order(coefs.biv)
labs <- c("Job Search","Contrib. to Entrep.","Garner Conts.","Personal Gain (% $100)","Group Gain (% $100)","Political Hom.",
          "Religious Hom.","Education Hom.","Class Hom.","Pers. Crisis","Pers. Success",
          "Prof. Crisis","Prof. Success","Least in Common","Most in Common","Pref. Interaction")

pdf("./1_Figures/SI_prot_causal.pdf")
  nf <- layout(matrix(c(1,2,3),1,3,byrow=T),widths=c(2,3,3))
  #layout.show(nf)
  par(mar=c(2,5,3,0))
  plot(rep(0,length(coefs.cause)), 1:length(coefs.cause), 
       type = "n", axes = F, xlab = "", ylab = "",main = "")
  axis(2, at = 1:length(coefs.cause), label = labs[ords], las = 1, tick = FALSE, cex.axis =.9,pos = .9)
  
  par(mar=c(2,0,3,1))
  pchs <- ifelse(abs(coefs.biv[ords]/ses.biv[ords]) < 2.33,19,21)
  cexs <- ifelse(abs(coefs.biv[ords]/ses.biv[ords]) < 2.33,.6,1.2)
  cols <- ifelse(abs(coefs.biv[ords]/ses.biv[ords]) < 2.33,rgb(0,0,0,.1),rgb(0,0,0,.8))
  ltys <- ifelse(abs(coefs.biv[ords]/ses.biv[ords]) < 2.33,1,1)
  bgs <- ifelse(abs(coefs.biv[ords]/ses.biv[ords]) < 2.33,NA,"white")
  plot(coefs.biv,1:length(coefs.biv),type='n',xlab="Coefficient",ylab="",yaxt='n',xlim=c(-.1,.2),main="Bivariate Regressions")
  segments(coefs.biv[ords] - 2*ses.biv[ords],1:length(coefs.biv),coefs.biv[ords] + 2*ses.biv[ords],col=cols,lty=ltys)
  points(coefs.biv[ords],1:length(coefs.biv),pch=pchs,col=cols,cex=cexs,bg=bgs)
  abline(v=0,lty=2)
  
  pchs <- ifelse(abs(coefs.cause[ords]/ses.cause[ords]) < 2.33,19,21)
  cexs <- ifelse(abs(coefs.cause[ords]/ses.cause[ords]) < 2.33,.6,1.2)
  cols <- ifelse(abs(coefs.cause[ords]/ses.cause[ords]) < 2.33,rgb(0,0,0,.1),rgb(0,0,0,.8))
  ltys <- ifelse(abs(coefs.cause[ords]/ses.cause[ords]) < 2.33,1,1)
  bgs <- ifelse(abs(coefs.cause[ords]/ses.cause[ords]) < 2.33,NA,"white")
  plot(coefs.cause,1:length(coefs.cause),type='n',xlab="Coefficient",ylab="",yaxt='n',xlim=c(-.15,.15),main="SEM Analysis")
  segments(coefs.cause[ords] - 2*ses.cause[ords],1:length(coefs.cause),coefs.cause[ords] + 2*ses.cause[ords],col=cols,lty=ltys)
  points(coefs.cause[ords],1:length(coefs.cause),pch=pchs,col=cols,cex=cexs,bg=bgs)
  abline(v=0,lty=2)
dev.off()
